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At a time when galaxy surveys and other observations are reaching unprecedented sky coverage 
and precision it seems timely to investigate the effects of general relativistic nonlinear dynamics on 
the growth of structures and on observations. Analytic inhomogeneous cosmological models are an 
indispensable way of investigating and understanding these effects in a simplified context. 

In this paper, we develop exact inhomogeneous solutions of general relativity with pressureless 
matter (dust, describing cold dark matter) and cosmological constant A, which can be used to model 
an arbitrary initial matter distribution along one line of sight. In particular, we consider the second 
class Szekeres models with A and split their dynamics into a flat ACDM background and exact 
nonlinear inhomogeneities, obtaining several new results. One single metric function Z describes 
the deviation from the background. We show that F, the time dependent part of Z, satisfies the 
familiar linear differential equation for 5, the first-order density perturbation of dust, with the usual 
growing and decaying modes. In the limit of small perturbations, S ~ F as expected, and the 
growth of inhomogeneities links up exactly with standard perturbation theory. In particular, we 
exhibit an exact conserved curvature variable, necessary for the existence of the growing mode, 
which is the nonlinear extension of the first-order curvature perturbation. We provide analytic 
expressions for the exact nonlinear 5 and the growth factor in our models. For the case of over- 
densities we find that, depending on the initial conditions, the growing mode may or may not lead 
to a pancake singularity, analogous to a Zel'dovich pancake. This is in contrast with the A = 
pure Einstein-de-Sitter background where, at any given point in comoving (Lagrangian) coordinates 
pancakes will always occur. Analyzing the covariant variables associated with the space-time, we 
derive the associated dynamical system, which we are able to decouple and reduce to two differential 
equations, one for J7a representing the background dynamics and one for S describing the dynamics 
of the inhomogeneities. Our models are Petrov type D, which we show by explicitly deriving the 
only nonzero Weyl scalar ^2, which does not depend on A. Since this is the only Weyl contribution 
to the geodesic deviation equation, A can only contribute to lensing through its contribution to the 
background expansion. 



I. INTRODUCTION 

From the very beginning of modern cosmology (see 
Einstein in [l| and the account in @), the guiding idea 
behind the construction of models has been what later 
became known as the Cosmological Principle (cf. 
and Rcfs. therein): the assumption that the Universe is, 
at any given time, homogeneous and isotropic on large 
scales is translated mathematically into a Robertson- 
Walker metric, i.e. a metric that is assumed to repre- 
sent a space average and it is therefore exactly homo- 
geneous and isotropic. In addition, the nontrivial hy- 
pothesis is made that this metric should be a solution of 
Einstein's equations, thereby giving rise to a Fricdmann- 
Lemaitre-Robcrtson- Walker (FLRW) universe model. In 
other words, general relativity (GR) is assumed to be 
the correct theory describing gravitational interactions 
between galaxies and clusters of galaxies. While it fol- 
lows logically that a Universe that is statistically homo- 
geneous and isotropic should be described on average by 
the Robertson- Walker metric, it is not at all obvious that 
this metric should satisfy Einstein's equations: given the 
nonlinear nature of the latter, averaging the equations 
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is not the same as considering the equations satisfied by 
the average (see |rll for this "averaging problem" ) . 

The Universe is however inhomogeneous. Even if it be- 
comes homogeneous above a certain scale (see e.g. fiof). 
the size of this "fair sample" 0] depends on how we mea- 
sure it and what we mean by "statistical homogeneity" 
(cf . [ll[ and [l2| for alternative views) . In any case, hav- 
ing assumed the Cosmological Principle, the growth of 
inhomogeneities and their effects are typically modeled 
with perturbation theory about a "background" FLRW 
model. Within this framework, the formation of non- 
linear structures at smaller scales is treated assuming a 
Newtonian approximation, with N-body simulations [13| . 
Most observations are interpreted assuming this Fried- 
mannian framework; in particular, distances are com- 
puted assuming a FLRW distance-redshift relation, i.e. 
completely neglecting inhomogeneities (cf. [l4[ for an at- 
tempt to include them). 

In the last three decades, the combination of cos- 
mic microwave background (CMB) radiation, large scale 
structure (LSS) and supernova type la observations has 
provided support for a flat FLRW universe model cur- 
rently undergoing an accelerated phase under the action 
of a dark energy c omp onent (l5l . Il6j . Consequently, the 
flat ACDM model [H G3 has emerged as the standard 
concordance model of cosmology [Iill2(| . GR is assumed 
to be the correct theory of gravity, with cold dark mat- 
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tcr (CDM) and a cosmological constant A dominating 
the dynamics of the late Universe. CDM is responsible 
for structure formation, with A making up the balance to 
have a spatially flat Universe and driving what is inferred 
to be, within this conceptual framework, cosmic accelera- 
tion. Currently no other model can explain observations 
equally well while being statistically favored, simply be- 
cause any alternative introduces more parameters, with 
not enough advantages (e.g. see [2l|-[24j]). 

Einstein's equations were constructed to automati- 
cally satisfy the energy-momentum conservation equa- 
tion, thus they must include a A term for full mathemat- 
ically generality Therefore, one can argue [25| that A 
is a free parameter of the theory to be determined by ob- 
servations, to be treated on the same footing as the cou- 
pling of gravity to matter - Newton's constant G. How- 
ever, this view is not very popular and many alternatives 
to A are considered in the literature @, HI, [TH, [2(| |27| in 
order to explain the effects that, within the homogeneity 
and isotropy assumption, are interpreted as acceleration 

Wc may divide these alternatives into two main groups, 
each with two subgroups. In the first, the Cosmological 
Principle is maintained, thus there is an acceleration but 
the cause is not a cosmological constant: i) either A is re- 
placed by an unknown dark energy fl6l [20 . \2l\ (or both 
CDM and A are replaced by unified dark matter, see 
[29M32| | and Refs. therein) or, ii) an alternative theory of 
gravity is assumed [H, [3Jj . The second alternative is to 
maintain GR and i) either consider possible dynamical 
back-reaction effects of structure formation on the over- 
all expansion, in order to construct an average model 
with acceleration [(| [jf. or ii) consider inhomogeneous 
exact solutions of GR [35|, [36( and study the effects of 
inhomogeneities on observations Q. In the latter case, 
the idea is that what is normally interpreted as accel- 
eration is at least partly due to observational effects of 
inhomogeneities. It turns out that, within very special 
models, one can even eliminate A entirely, for instance 
for observers near the center of a spherically symmetric 
model (see [37j and references therein). However, it is an 
open question whether this can be achieved in inhomo- 
geneous models which exhibit an average homogeneity 
on the scales on which it seems to be observed [lOj, i.e. 
models in which the cosmological principle is preserved. 

Given the above arguments, it is of great interest to 
consider simple models where nonlinearities and inho- 
mogeneities are fully taken into account while, at the 
same time, a FLRW background can be clearly identified. 
A very interesting class of models was introduced long 
ago by Lindquist and Wheeler [38j and revised by Red- 
mount [1^], which has then been recently reconsidered 
and generalized by Clifton and Ferreira |40( , where point- 
like masses (represented by Schwarzschild black holes) 
are distributed in a lattice and the overall expansion is 
described by the Fricdmann equation. The main moti- 
vation of Clifton and Ferreira in following the Lindquist 
and Wheeler construction is the observation that indeed 



the Universe largely consists of galaxies and clusters of 
galaxies surrounded by vacuum. The question they ad- 
dress is how observations and measurements of the cos- 
mological parameters are affected in a highly inhomoge- 
neous universe whose overall dynamics is homogeneous 
and isotropic. However, this lattice construction is only 
an approximate solution to Einstein's equations, with "no 
man's land" in between the matched spheres and, per- 
haps even more importantly, where the inhomogeneities 
are strongly nonlinear at all times. 

Having in mind the same type of questions addressed 
in (4^1 , in this paper we adopt a milder approach to the 
construction of models of nonlinear inhomogeneities in 
a clearly identifiable FLRW background. Specifically, we 
consider exact solutions of GR where, starting from stan- 
dard small perturbations of a FLRW universe, the mat- 
ter distribution is continuous and can evolve to a highly 
nonlinear stage. In the process, the inhomogeneities can 
either form a distribution of large voids or over-densities, 
or a mixture of the two, with over-densities possibly even 
forming pancakes as in the Zcl'dovich approximation in 
Newtonian cosmology. The benefit of our model is there- 
fore two-fold: i) we consider exact solutions of Einstein's 
equations, therefore avoiding any possible problem asso- 
ciated with approximations and matching and ii) these 
exact solutions describe nonlinear inhomogeneities grow- 
ing on top of a FLRW background, with the possibility 
of modeling a rather arbitrary distribution of both voids 
and over-densities. It should be noted that there is no 
restriction on the inhomogeneities to have a zero average 
in any sense, thus the FLRW background may or may not 
be representative of an average; we shall not particularly 
address this issue here, leaving it for a future analysis 

a 

We do not necessarily expect that considering observa- 
tional effects of nonlinear inhomogeneities can eliminate 
- entirely or in part - the need for dark energy. On the 
contrary, the study of the nonlinear growth of perturba- 
tions and the effects they may have on observations is 
a topic of the greatest interest in the context of ACDM 
cosmology. Therefore we include a cosmological constant 
A in our models. The exact solutions we are considering 
are a generalization with A of the pure dust models of 
Szekeres [42j and arc perse not new, as they were first 
found by Barrow and Stein-Schabes [43| , who used them 
to prove a nonlinear version of the cosmic no-hair the- 
orcixQ. However, what is new is the form in which we 
present this solution, which allows for a clear distinction 
between inhomogeneities and the FLRW background. In 
doing this, we generalize to the A / case the work 



1 The cosmic no-hair conjecture states that, for all solutions of 
Einstein's equations with A > and matter satisfying the weak 
energy condition, the future asymptotic state is stationary inside 
the cosmological event horizon of any future inextendible timelike 
curve. Dynamically, this states that de Sitter space-time is the 
future attractor; see [43, 44] and references therein. 
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by Goode and Wainwright (45[, focusing on the particu- 
lar subclass of solutions that admit a flat FLRW back- 
ground. Through the split between FLRW background 
and inhomogeneities, we achieve many new results. In 
this paper, we present the solutions and their properties, 
while we leave the analysis of the observational effects of 
inhomogeneities in these models to a second paper [4lj |. 

In order to avoid any confusion with approximate so- 
lutions, in this paper we refer to deviations (from the 
background) rather than perturbations, even in the case 
of standard variables such as the dimensionless density 
inhomogcneity 6 = (p — p)/p, where p and p respectively 
indicate the density in the inhomogeneous universe and 
in the FLRW background. We choose the latter to be 
spatially flat, see Eq. ([15]) for the metric where, in syn- 
chronous comoving coordinates, Z is the single metric 
function that characterizes the deviations from this back- 
ground. The main results of our analysis are: 

• generalizing the Goode and Wainwright [45[ formal- 
ism to the A ^ cas^l, we exhibit exact solutions 
for the growing and decaying modes of the metric 
deviation Z, assuming a flat ACDM background; 

• as for the A = case, the second order ordinary 
differential equation (ODE) for the growing and de- 
caying modes of Z is linear, and is the same equa- 
tion that it is satisfied by 5 in the linear regime; 

• therefore, as for 8 in the linear regime in a flat 
ACDM background, the growing mode in Z asymp- 
totically approaches a constant value; 

• given that the equation for the metric deviation Z 
is linear, Z satisfies a superposition principle even 
in the nonlinear regime, while 5 does not; 

• we explicitly show that the second order ODE sat- 
isfied by Z admits a first integral; this results in a 
conserved curvature deviation, d, which turns out 
to be strictly related to the growing mode, i.e. d = 
implies a purely decaying mode; 

• we eliminate the residual gauge freedom from the 
form of the metric in [45j |. so that Z = 1 immedi- 
ately corresponds to the FLRW background; 

• we exhibit the exact analytic solution for the non- 
linear 5, showing the analogy to the Newtonian 
Zel'dovich pancake for the case of over-densities; 

• we discuss the formation of singularities at any 
given comoving point0 and find that, for the case 



2 In this paper we restrict our generalization to the second class of 
Szekeres-like metrics considered in [45| . 

3 Our synchronous comoving coordinates should be thought of as 
the equivalent of Lagrangian coordinates in a Newtonian descrip- 
tion. 



of over-densities, pancakes may form in the future, 
while the past singularity is a FLRW-like Big-Bang 
for purely growing modes in Z, and spindlclikc 
(cigar, or Kasner-like) if a decaying mode is present 
in Z; 

• given that the growing mode in Z evolves to a con- 
stant, at any given point pancakes do not form if 
appropriate initial conditions are chosen at that 
point, unlike the A = case where pancakes are 
unavoidable; 

• we explicitly work out the analogy between our for- 
malism and that of perturbation theory, exploiting 
the linearity of the equation for Z and the con- 
served curvature variable d; the latter plays the role 
of the conserved linear curvature perturbation; 

• we explicitly give an analytic expression for the 
growth factor of the density deviation 5, valid into 
the nonlinear regime; 

• we reconsider the d yna mical system characterizing 
our models (cf. j46H48| ) in terms of density p, ex- 
pansion O, shear a a b and Electric Weyl tensor E a b 
and, introducing new dimensionless variables, we 
decouple the system of ODEs, reducing the dy- 
namics to 2 equations, one for describing the 
background and one for 6 describing the inhomo- 
geneities; 

• we verify that our solutions are of Petrov type D 
[49j , explicitly finding the null tetrad that makes ^2 
the only nonvanishing Weyl scalar, which we show 
does not directly depend on A: it follows that in our 
models A can only affect lensing indirectly, through 
the coupling of the dynamics of the background to 
that of the inhomogeneities. 

From this list, it is clear that some of our results are 
directly related to astrophysical cosmology, while others 
are more in the domain of rclativistic cosmology. The 
reader more interested in the former can first read Sees. 
ITT1 IIII Al and IIV| leaving the other sections for a second 
reading. 

The plan of the paper is as follows. In Sec. [TT| we 
present a summary of our analysis of Einstein's field 
equations (EFEs); full details are given in Appendix 151 
After choosing a fiat ACDM background, we present ex- 
plicit solutions for the growing and decaying modes of 
the metric deviation Z . 

We give a physical interpretation to our models in Sec. 
IIII Al looking at what kind of nonlinear density distribu- 
tions are possible. In Sec. IIII B I we present an analysis of 
the singularities: in particular, we look at when pancakes 
form or do not form, depending on the initial conditions. 
The space-time is analyzed using the Petrov classification 
scheme in Sec. IIII CI 

In Sec. IIVI we analyze the analogy between our exact 
solutions and cosmological perturbation theory. In par- 
ticular, in our model density inhomogeneities can grow 
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highly nonlinear; however, we would like to use initial 
conditions as they are used in standard perturbation the- 
ory and hence we show how, in the linear regime, the two 
are related. 

In Sec. El we consider the dynamical system associated 
with a set of covariant variables. First, in Sec. IV Al we 
relate the covariant variables with the metric functions. 
Then, in Sec. IV Bl we show how the dynamical system 
for the covariant variables can be reduced to only two 
ODEs, and we present a phase plane analysis. 

In Appendix [A] we demonstrate how the continuity 
equation directly implies the form of the density and den- 
sity deviation in our model. In Appendix [Bl we give de- 
tails on how we solved the EFEs, with the dimensions of 
all the variables and parameters given in Appendix [D] 
Finally, in Appendix [El we demonstrate how we can ob- 
tain axial symmetry in our model, depending on how we 
choose the free functions. 

Throughout the paper, we choose units c = 8irG = 1. 



II. SOLVING THE EFES 
A. Setup 

In this paper we shall consider the second class 
Szekeres-type metrics 42j, in the form that Goode and 
Wainwright H^] introduced 



It should already be clear that S = S(t) in the metric 
([1]) acquires the role of a FLRW scale factor. In any case, 
a simple interpretation of the metric ([I]) is immediately 
obtained if we consider the generalization of the Hubble 
expansion law 0, [5^] ■ Consider two fundamental comov- 
ing observers moving with 4-velocity u a and connected, 
at any given time t, by a vector X a (thus orthogonal to 
u a ). For the components of X a we find 



HX 2 



X z = \H + - \ X z 



(6a) 
(6b) 



where H = S/S. Hence we deduce that the Hubble law 
along the x- and y-axis is unmodified, i.e. the same as 
that of a FLRW background, whereas along the z-axis it 
is changed by the inhomogeneities encoded in Z . As we 
shall see in Section [Vl 6 = Z/Z is precisely the deviation 
of the expansion (scalar) from that of the background, H. 
Thus in our synchronous comoving representation, which 
is the relativistic analog of a Newtonian Lagrangian de- 
scription, fluid elements occupy a fixed coordinate posi- 
tion, but the physical distance between any pair of them 
along the z-axis is modified by Z. As we shall see, this 
may even lead to pancakes, when Z — > 0, in analogy 
of the Zel'dovich pancakes of the Newtonian-Lagrangian 
description. 



ds 2 = -dt 2 + S 2 



e 2a(x) ^ x 2 + dy 2j + Z 2 dz , 



(1) 



where S = S(t) and Z = Z(x,t). We will generalize 
their analysis, including a cosmological constant A in Ein- 
stein's equations [43j : 



G a b — T a b — h-g a b- 



(2) 



The matter content, a pressureless dust component, rep- 
resents CDM. Fluid elements move along geodesies and, 
in the synchronous coordinates in the metric ([1]), these 
geodesic flow lines are orthogonal to the cosmic time t 
hyper-surfaces, with 4-velocity 



[1,0,0,0]. 



(3) 



The coordinates in Eq. ([I]) are therefore also comoving 
and the fluid flow is irrotational. The energy-momentum 
tensor Tab only has one nonzero component, which is 



rpU 



(4) 



where p is the energy-density of the dust. It follows di- 
rectly from the conservation equation T ab '-b = (see Ap- 
pendix [3J that in general for this metric 



M(x) 
S 3 Z ' 



(5) 



where for now M is a general function of space, which 
will be restricted by the EFEs later on. 



B. Summary of the calculations 

The details of our analysis of the EFEs for the metric 
(fT|) are given in Appendix [Bl From this analysis, we 
obtain the following equation for the dimcnsionless scale 
factoB S{t): 



3 P 



1 A K 

3 A+ ^' 



(7) 



which we recognize as the Fricdmann constraint equa- 
tion for ACDM, where we have defined the homogeneous 
energy-density 



S 3 ' 



(8) 



Po and Sq = 1 are the values of p and S today. The curva- 
ture constant K has dimensions and is either vanishing, 



4 In the main part of the paper we use an over-dot to denote dif- 
ferentiation with respect to t. In the Appendix, however, we use 
a subscript t, for uniformity with the space derivatives in EFEs. 
The latter are completely solved with respect to space coordi- 
nates, hence only ODEs with respect to t need to be analyzed in 
the main part of the paper. 
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positive or negative for a flat, closed or open universe re- 
spectively (see e.g. Q); it is linked to the metric through 
the relation 



(9) 



The function Z in the line element ([T]) can be split as 

Z(x,t)=F(z,t) + A(x), (fO) 
where A can be written in the forrrH 

A(x) = a(z) + b(z)x + c(z)y + d(z)(x 2 + y 2 ), (11) 
and F obeys the following linear homogeneous ODE 



f + 2|f-|f = o. 



(12) 



This equation is well known: as it was noted in |45| . it 
is indeed the equation satisfied by the first-order density 
perturbation in a dust (CDM with or without A) FLRW 
universe (see for instance [2|). Less well known is that 
this equation admits a first integral (see Appendix [C] for 
a derivation) 



S ■ p 2d 
—F +-F =0 

5 2 S 2 ' 



(13) 



where, through the field equations, the conserved quan- 
tity turns out to be the curvature variable d, appearing 
in A, Eq. (fTT|) . Clearly Eq. (fT2|) admits two linearly in- 
dependent solutions: 



F(z,t)=0+(z)f+{t) + P-(z)f-(t), 



(14) 



where f+(t) represents the so called growing mode and 
f—(t) the decaying mode. Alternatively, f-(t) is the so- 
lution to the homogeneous part of Eq. (fl"3"j) and f+(i) 
is the particular solution. We will explore in more de- 
tail the significance of the functions 5, /+, /_ and the 
conserved quantity d in the following sections. 



C. Solving for the metric functions 

We focus our attention on the K = case (i.e. a flat 
background) and therefore, the metric can be written in 
the form 

ds 2 = -dt 2 + S(tf [dx 2 + dy 2 + Z(x, tfdz 2 } . (15) 
Using Eq. J8|) , Eq. (jTj) for the K = case gives 



q2 _ P0 . A 2 

b ~ 35 + 3 b ' 



(16) 



5 Note that this expression for A will not be used in following 
Sections, since it contains gauge functions, which will be fixed in 
SeclHCl 



which enables us to embed our space-time (fTB"]) in a 
FLRW background through the scale factor S(t). There- 
fore the metric (fT5"j) can be seen as describing an exact 
inhomogeneity, specified by Z, in the ACDM background 
described by S(t) and parameterized by fl m = po/(3H 2 ) 
and SI a — A/(3Hq), where Hq is the Hubble parameter 
and fl m = 1 — f^A- Solving Eq. (|16p then gives 



S(t) = 



1 - 0/ 



1/3 



sinh 2/3 ( -HoVSht), 



where we have set to zero an integration constant which 
would only shift the time of the initial singularity. With 
this, the age of the Universe today is 



in 



2H n 



: arcoth fi 



.-1/2 
'A 



(18) 



it is easy to check that inserting this expression into Eq. 
([171) yields 5(0) = 5 = 1, as it should. Now, defin- 
ing the dimensionless variable r = 
equation ([T2")) for F simplifies to 



^t, the differential 



F" + ^coth(r)F'-^— 



3 sinh (r 



■F = 0, 



(19) 



where a dash denotes the derivative with respect to r. 
The two linearly independent solutions are 



/- = 



cosh(r) 
sinh(r) 



cosh(r) f sinh 2/3 (r) , 
./+ = . . , ; / . o . . dr. 



sinh(r) J cosh (r) 



(20a) 
(20b) 



The function S(t) in Eq. (|TT|) is the scale factor of 
the ACDM background model. The differential equation 
(fT9| . obeyed by F, is exactly the one obeyed by the first- 
order density perturbation in this background (for an in 
depth discussion, see e.g. 0). The two independent so- 
lutions for F are shown in Fig. (fTJ). 

It should be noted that the same formalism applies in 
the A = case [45|, where the FLRW background is an 
Einstcin-de-Sitter model. The main difference is that the 
growing mode /+ grows linearly with the scale factor 5 
for A = 0, whereas in our case (see top panel in Fig. (Q])), 
/+ asymptotically approaches a constant value. 

Now, we would like to remove some residual arbitrari- 
ness in the free functions of z appearing in Z (a, b, c, d, 
j3 + and First, Eq. (fL3|) can be used to express the 

conserved curvature d as a function of F and 5: 



55 • 

— F + —F. 

2 45 



Pa 



(21) 



Then, one can substitute the general solution for F, 
which has been computed earlier; after some algebra one 
finds 



d(z) = B0+(z), 



(22) 
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FIG. 1. Plots of the growing (top panel) and decaying (bot- 
tom panel) modes of the solution for F, as derived from Eq. 
(|12|l . The solutions plotted here are given in Eqs. (I20a|l and 
(l20bl). 



where we have defined 

B= 1 -(p 2 A) 1/3 = 3 -H 2 [n A (i-n A ) 2 ] 



1/3 



(23) 



as the constant which gives dimensions of L 2 to d. This 
can be substituted into A to obtain 

A{x) = a(z) + b(z)x + c(z)y + B[3 + {z)(x 2 + y 2 ). (24) 

Hence there are still 5 free functions of z remaining in 
the line element (a, b, c, /3 + and /?_). This freedom can 
be further reduced by one order by transforming the z- 
coordinate. First of all, however, let us introduce two 
new functions 7 and uj; with no loss of generality, we can 
write 



b(z) = 2^(z)Bj3 + (z), 
c(z) = 2u(z)Bp+(z). 



(25a) 
(25b) 

Then, substituting these in Eq. (|24|) we find 

A = a + B/3+ [{x + 7 ) 2 + (y + uj) 2 - ( 7 2 + w 2 )] . (26) 

Now we have to make a transformation in the z coordi- 
nate and rescale /3+ and /3_. We choose 



/3_(z) 



z = I {a(z) - B/3 + (z) [ 7 2 (z) + u 2 (z)]}dz, (27) 



a(z)-B(3 + (z) [^(z)+^(z)Y 



a(z)-B0 + (z) [ 1 2 (z)+o J 2 (z)Y 



(28) 
(29) 



With this coordinate transformation and rescaling of /3+ 
and (3- , we obtain the simplification 



\x + l{z)f + [y + oj(S)} 2 } 



(30) 



A(x,y,z) = l + Bp + (z){[ 
and 

Z(t, x, y, z) = P+(z)f+(t) + /L (2)/_ (t) + A{z). (31) 

With these we can now drop all the tildes and write our 
metric in the final form 



ds 2 



-dt 2 + S{t) 2 



rf.r 



+ dy 2 + (l + P+(z)f+(t) + /8_ («)/_(*) + Bp + {z) [[x + 7 (^)] 2 + [y + u(z)} 2 }) 2 dz' 



(32) 



We have thus reduced the freedom in the metric to four 
free functions (7, us, j3+ and /?_). This expression also 
clarifies the meaning of the coordinate transformation we 
have just performed: in essence we have completely fixed 
the gauge, such that when /3+ and /3_ are equal to zero, 
which implies F = and Z = 1, our metric exactly takes 
the form of the background FLRW space-time. Given 



I 

these arguments, we will henceforth be using Eq. (|30p as 
the expression for A(x, y, z). 
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III. INTERPRETING AND CLASSIFYING THE 
SOLUTION 



A. Visualization of simple solutions 



In this section we would like to give some intuitive un- 
derstanding of what kind of energy-density distributions 
are possible in the developed space-time. We can see 
from Eq. (|aT0]) and ((B28|) that 



A(x) 



M(x) 



Pi) 



(33) 



therefore we can rewrite expression ((5]) for the density, 
eliminating the dependence on M 



P 



AI 



S 3 Z S 3 (F + A)' 



(34) 



Defining as usual the dimensionless density deviation 
from the background density p (Eq. ([8])) 



6 = 



we obtain 



F 



A 



F 
Z' 



(35) 



(36) 



Using the decompositions of F and A derived earlier, we 
can write 



5(x,y, z,t) 



1 + P+{z)f+(t) + £-(*)/-(*) + BP+(z) {[x + 7 (z)] 2 + [y + Lu{z)f} 



(37) 



In general, the density field in (|37| can represent, at any 
given space point x, either an over-density or an under- 
density, depending on the values of f3+ and /3_. However, 
for the case of over-densities, at any given time, there 
exist points where (|37p necessarily implies a pancake sin- 
gularity. In essence, this is due to the vanishing of the 
function Z in the denominator of (|36[) - (f3T|) . We will dis- 
cuss the existence and properties of these singularities in 
more detail in the next section. However, it is important 
to note here already that these pancake singularities are 
only due to the continuous description of matter in our 
models, i.e. they are shell crossing singularities (see e.g. 
[HI) and are analogous to Zel'dovich pancakes in New- 
tonian gravity. 

The structure of the density distribution (|37[) is hard 
to visualize in the general case 7 ^ and/or lu 7^ 0, 
which also implies the absence of Killing vectors for our 
space-time. Purely for illustrative purposes, we now con- 
sider the restricted case 7 = lu = 0, which implies axial 
symmetry, see Appendix [El 

As a first example of the inhomogeneous density distri- 
bution (j37|). we consider in Fig. ((2]) a purely under-dense 
growing mode, /3 + > 0, with /3_ = 7 = lu = 0. In this 
case there are no pancakes and the density field is reg- 
ular everywhere. The function j3+ has been chosen to 
take the form j3 + cx [1 — sin(fcz)] for k = IMpc -1 and we 
are plotting the density deviations at an arbitrary value 
of t. On this plot, one can see that the center of inho- 
mogeneities runs along the z-axis, which is due to the 
condition 7 = lu = 0. In general, this center can take any 
path around the z-axis, set by the two functions 7(2) and 

LU{Z). 



A distribution of over- and under-densities (or pure 
over-densities) can be simulated, if one only considers 
the space-time "close" to the z-axis, in order to avoid 
pancakes. An example of this is shown in Fig. <j3j> . We 
will discuss this case in more detail in the next section. 




r[Mpc] ° 10 z[Mpc] 

FIG. 2. The 8 profile of an under-density at an arbitrary 
time ( in a ACDM background with J1a = 0.75 and Ho = 
72kms- 1 Mpc- 1 . Wc assume a purely growing mode with 
P+ oc [1 - sin(fcz)] for k = lMpc -1 and /3_ = 7 = uj = 0. 
In this case the space-time is axi ally sym metric, so r is the 
distance from the z-axis, r = \/ x 1 + y 2 . All distances are 
comoving and given in Mpc. 
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r[Mpc] 10 z[Mpc] 

FIG. 3. The S profile for a distribution of over- and under- 
densities at an arbitrary time t in a ACDM background with 
= 0.75 and Hq = 72 kms _1 Mpc _1 . We assume a purely 
growing mode with /3+ oc sin(fcz) for k = IMpc -1 and /3_ = 
7 = uj = 0. As in Fig. J2j, r is the distance from the z- 
axis, r = \/ x 2 + y 2 . The increasing behavior of the over- 
densities away from the z-axis is due to singularities located at 
a certain r*, beyond the boundaries of the plot. All distances 
are comoving and given in Mpc. 

B. Classification of singularities 

The cosmological model we are considering contains 
only irrotational dust: singularities in these space-times 
have been studied for a long time HU, see also [H)]. 
Therefore, the types of singularities we encounter in our 
model are well known, however it is essential to under- 
stand if and when they occur. 

Considering the fact that /+, /_ > 0, we can see from 
Eq. (f3"Tl) . that if 0+(z) and 0-(z) assume some negative 
values, we have a singularity in 8 if the denominator (Z) 
goes to zero. This also causes a singularity in p. This is- 
sue has been considered in great detail for the case with- 
out a cosmological constant by Goode and Wainwright 
[52l ]. We define our metric functions slightly differently, 
such that 0+(z) and 0-(z) have opposite signs to the def- 
initions in [5 2( | . The analysis of the singularities for the 
case A 7^ gives different results than for A = 0, which 
can be exploited to model inhomogcncities in a physically 
meaningful way. 

First of all, let us introduce some formalisms. We de- 
fine the variables l a through 

3 

ds 2 = -dt 2 + ^2l 2 a (dx a )\ (38) 

a=l 

where in our space-time we find that l\ = I2 = S and I3 = 
SZ. Using these variables, one can classify singularities 
into three different types (see e.g. |4^|): 

1. a pointlike singularity when all three l a — >• , 



2. a cigar or spindle singularity if two l a — > and the 
other one diverges, 

3. a pancake singularity if one l a — > and the other 
two converge to a finite value, 

as we approach the singularity. We remind the reader 
that we use synchronous comoving coordinates, so that 
in the following "fixed space point" refers to these coor- 
dinates. 

We find that the initial singularity at t = can only 
be either a pointlike singularity or a cigar singularity de- 
pending only on /?_ : if /?_ 7^ we get a cigar singularity 
and if 0— = we get a pointlike singularity. In other 
words, if a decaying mode is present then the initial sin- 
gularity is velocity dominated (see [Hlj], cf. also [53| and 
references therein) and Kasner-like, while when we only 
have a growing mode, the initial singularity is matter 
dominated and effectively isotropic: in approaching the 
singularity, the growing mode decreases and our space- 
time becomes FLRW with a small perturbation. 

As we argued earlier, there can also be singularities 
at some time > 0, where we find divergences in the 
density field. We find that these singularities can only 
be pancake singularities. We present the singularities in 
Table U for the case of 0- < and in Table HT1 for the case 
of /?_ > 0. 

To analyze the case 0- < 0, we need a more in depth 
understanding of the behavior of several functions. The 
two time dependent functions /_ and /+ exhibit an 
asymptotic behavior for large values of t, as it can be 
seen from Fig. [TJ We find /_ — > 1 and /+ — > /?° as 
t — > 00, where /?° is a finite and positive number. This 
is significantly different from the A — case, where /+ 
does not asymptote to a finite value for large t. Since /_ 
and /+ have an asymptotic value, it follows that for every 
fixed space point F and hence Z have a finite asymptotic 
value as well, see Eqs. (fTU)) and (|30[) for the definition of 
Z . We therefore introduce the new parameter Zoo, which 
is the asymptotic value of Z at a given space point. We 
would now like to deduce that Z always has a maximum 
value for /?_ < and 0+ < 0. At early enough times, the 
decaying mode /_ dominates over / + and so Z w /?_/_, 
which is positive, since 0- was assumed to be negative 
and /_ is a strictly decaying function. At late enough 
times, we can rearrange Eq. (fT3")) to give Z « 2B0 + /(SS), 
which is negative, since B, S and S are strictly positive 
and we have assumed 0+ to be negative. We therefore 
find that Z changes sign, from positive to negative and 
hence must have a maximum, which we call Zm- Ana- 
lyzing the parameters Zm and Z^ aids the distinction 
between different cases in Table HI 

In the case of 0- > we find a splitting between dif- 
ferent cases, depending on a new parameter. Since Z is 
positive initially (1 + /?-/- being positive and dominat- 
ing at early times), we find that if Z m is positive, Z has 
no zeros and hence we find no pancakes. This distinction 
between cases turns out to depend on the value of 0+. 
Clearly, if + is positive (and hence Z is positive for all t), 
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TABLE I. Classification of singularities occurring at some fi- 
nite time t„ > for /3_ < 0. 
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Zm 


Zoo 


number of pancake 
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< 
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2 


> o 




< 





> o 




> 


1 



we only have under-densities and hence no singularities. 
If /3+ is negative, we find a critical value, which divides 
the cases of singularities and no singularities. This value 
occurs when 

Zoo = l + |/3+| {/+ +5 [(x + 7 ) 2 + {y + lo) 2 } } , 

(39) 

vanishes. For a fixed space point (and hence given values 
of /3_, x, 7, y and to), we can find this critical value |/3 +| 
to be 



1/3- 



1 



f? + B[(X + 7 )2 + (y + U )i 



(40) 



This parameter is used in Table |TT] to decide whether a 
given point in space will have a future singularity. 



TABLE II. Classification of singularities occurring at some 
finite time U > for /3_ > 0. 



p+ 


number of pancake singularities 


> o 





<0, < 





<0, |/3+| > |/3+| 


1 



The second case in Table [TT1 where /3+ < and 
|/3+| < |/3+|, is the most interesting from a pragmatic 
point of view. Given a certain point in space, we can 
always find a value |/3+|, such that for |/3+| < |/3 + | there 
will be no future singularity at this point. This is dis- 
tinctly different from the A = case where this is not 
possible. In practice this means that if we would like to 
model some density distribution on a certain space re- 
gion without incurring in a pancake, we need to find the 
maximum |/3+| within this region, which will restrict the 
maximum over-density we can model. " Close" to the z- 
axis this restriction (for /3_ = 7 = uj) turns out to be 
fairly weak. For instance, using initial conditions at re- 
combination, we can start with initial values even greater 
than the measured power spectrum amplitude [28[ and 
have no future pancakes in a finite region around the 
z-axis. 



C. Petrov classification 

The Petrov classification is used to distinguish different 
types of space-time metrics by analyzing algebraic prop- 
erties of the Weyl tensor (for a discussion of the Weyl 
scalars and the Petrov classification, see [HI, HH). A 
main point to be noted is that these properties are purely 
geometrical and unrelated to the theory of gravity con- 
sidered. However, an understanding of these properties 
helps the physical interpretation, especially in those cases 
where the space-time can be seen as a nonlinear pertur- 
bation of some background. Such physical interpretation 
was given for instance by Szekeres in [56j . and is based 
on the so called Weyl scalars. These are defined by 



*o = C abcd k a m b k c m d , 
* 1 = C abcd k a m b k c l d , 
* 2 = -C abcd k a m b l c fh d 
*3 = C abcd l a fh b l c k d , 
*4 = C abcd l a m b l c m d , 



(41a) 
(41b) 
(41c) 
(41d) 
(41e) 



where C abc d is the Weyl tensor and k,l,m and m is a 
complex null tetrad. These five complex scalars repre- 
sent, in four dimensions, the ten components of the Weyl 
tensor. In essence, the Petrov classification of a certain 
space-time involves finding the complex null tetrad such 
that the number of Weyl scalars reduce to a minimal set. 
If one can find a tetrad such that the only non- vanishing 
Weyl scalar is <3/2> then the space-time is said to be 
Petrov type D. Well known examples arc Schwarzschild 
and Kerr. A space-time containing gravitational waves 
necessarily has nonzero \&o an d ^4, e.g. a perturbed Kerr 
[H, [57], HH]. For this reason, an analysis of the Weyl 
scalars can be used in extracting gravitational waves in 



On the other hand, 



numerical relativity, see e.g. 
^0 and \&4 can be nonzero in a space-time with no gravi- 
tational waves, e.g. as in the case of a stationary rotating 
star [62|]. Clearly, a Petrov type D space-time does not 
contain any gravitational radiation. 

The general result that Szekeres space-times are Petrov 
type D is well known, see e.g. (43, [Hf, therefore there 
are no gravitational waves. Naively, without the knowl- 
edge of the Petrov type and its meaning, this is counter- 
intuitive: since we have time dependent matter inhomo- 
geneities, one would expect gravitational radiation to be 
present. Since our metric (|T5l) has the Szekeres form, it 
must also be Petrov type D. We now want to show this 
explicitly, especially to analyze the FLKW limit of our 
model. For a derivation of the complex null tetrad, see 
Section IV Al Using this basis, we obtain that the only 



6 Another way to conclude that no gravitational radiation is 
present in Szekeres space-times is that the the magnetic Weyl 



tensor vanishes 1491 1. see Section |yj 
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nonzero Weyl scalar itQ 



$0 




(42) 



This expression is derived from the metric alone and so 
contains only geometric information. In particular this 
shows that our space-time has a single independent Weyl 
component. Using the EFEs, we can now relate this ex- 
pression to the matter content and we find 



(43) 



2 6S 3 Z 6S 3 6 P 



In the case of under-densities, it follows from Eq. (|37|) 
that this expression goes to zero for large x 2 + y 2 and 
large t. In the case of over-densities, at any given space 
point, ^2 w iU diverge when a pancake forms, except in 
the second case in Table [III where there is no pancake 
and \I/2 — >• for large t. However, in approaching the 
pancake at t* it turns out that 



_*2 

92 



lpo-FZ. 
6 S 3 Z 2 ' 



(44) 



see Eq. (|57j) below for a definition of O. At t* all quan- 
tities in this expression have a finite value, while Z = 0; 
therefore this dimcnsionlcss measure of the Weyl curva- 
ture vanishes at the pancake. A space-time with only a 
nonzero is Petrov type D and a space-time with all 
Weyl scalars identically zero is type O, i.e. conformally 
flat. This means that we have a type D space-time in 
general. For all cases without a pancake, the space-time 
tends to a type O and FLRW space-time for large values 
of x 2 + y 2 or large values of t. We notice from Eq. (|43|) 
that ^2 does not explicitly contain A. Since ^2 is the 
only Weyl contribution to the geodesic deviation equa- 
tion (60l |. this shows that there is no direct contribution 
to lensing from the cosmological constant through the 
Weyl curvatur^]. 



IV. RELATION TO PERTURBATION THEORY 

Cosmological perturbation theory concerns itself with 
the dynamics of small deviations from a homogeneous 
FLRW background and the corresponding approximate 
treatment of Einstein's equations. In the covariant ap- 
proach to perturbation theory [64| . the variable A is in- 
troduced to analyze the behavior of density perturba- 
tions, see [IHHH. This gauge-invariant variable reduces 
to the density perturbation 5 in the comoving gauge and 



7 Z X x represents the second derivative of Z with respect to the x 
coordinate, see Appendix [B] 

8 The issue of a direct contribution from A to gravitational lensing 
has been the subject of recent investigations, see e.g. [631] and 
references therein. 



to the corresponding gauge-invariant variable A derived, 
in Fourier space and within the metric perturbation ap- 
proach, by BardeerU [§3] and Kodama-Sasaki (68[. It can 
also be used in perturbative studies of the behavior of 
inhomogeneities and anisotropies in the neighborhood of 
isotropic singularities for more general space-times such 
as Bianchi models, see e.g. [69L r70l]. 

The second order differential equation governing the 
evolution of A for prcssureless dust is 



A + 2§A-ip-A 



0. 



(45) 



The same equation is satisfied by S in Newtonian theory 
0. In general, the dynamical content of the second order 
equation (|4"5"|) can be reexpressed by a system of two first- 
order equations, coupling A to either CorZ, due to the 
constraint 



C = -4SSZ + 2S 2 pA, 



(46) 



where Z and C represent the spatial variation of the ex- 
pansion scalar (|57| and the 3-Ricci scalar respectively 
[651 ] . Since we are considering irrotational dust, the sys- 
tem for A and C takes the form 



A 

C 



c_ 

S 2 ' 

0, 



(47) 



therefore, clearly, C is a quantity representing a con- 
served curvature perturbation, with dimensions L~ 2 . 

For the metric developed in this paper, we find the 
background expansion to be given by 5". The other time 
dependent function F represents the deviation from ho- 
mogeneity. The second order differential equation which 
describes its evolution is derived in Appcndix[B]and given 
in Eq. (fT2j) . We now notice that this is exactly the same 
equation as the one satisfied by A in Eq. (|4"5")) . The same 
tactic of reducing the second order equation to a set of 
two first-order equations can then be employed for this 
variable and we find 



F 



2B/3+ 



S 2 ' 

(2Bf3 + Y=0, 



(48) 



where B is the constant defined in Eq. ([23]) . From a 
dynamical system perspective, it is interesting to note 
that the first equations in the systems (|47p and ([4"5[) arc, 
respectively, first integrals of Eq. (j45[) and Eq. (|T2|) . with 
C and 2Bf3 + the corresponding conserved quantities; see 
Appendix [C] for the explicit integration. 

It is striking that the differential equations for A and F 
take the same form. We can make the analogy even more 
apparent by considering the limit in which 8 is small. 
From Eq. ([36]) we can solve for F in general 



F = — 



AS 

S + V 



(49) 



In Bardccn notation this variable IS dri 
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so that for small values of S, F « —AS. This can be 
substituted into the differential equation for F to obtain 



Q 



3*' 

= 0, 



(50) 



i.e. system f|47[) . where we have defined Q — — 2 ^f + ■ In 
the limit of small S we have thus retrieved the differential 
equations governing the growth of density perturbations 
in cosmological perturbation theory, as one would expect. 
However, the strength of our model is that Eqs. (fT2|) or 
(I48[l can be used to evolve F into the nonlinear regime, 
with Eq. (|36p giving the corresponding S. Moreover, we 
have not only found the analogy of the growth of per- 
turbations to cosmological perturbation theory, but we 
have also found a conserved quantity into the nonlinear 
regime in our system - 2Bf3 + . The systems of differen- 
tial equations admit a decaying solution only for a zero 
conserved quantity, C and 2B(3 + respectively for the lin- 
ear and exact nonlinear regimes. This shows that the 
growing mode solution for F or S corresponds to a par- 
ticular solution to the respective equations, generated by 
a nonzero conserved curvature inhomogencity, either the 
exact 2Bf3 + or the perturbation C, respectively. 

Finally, using (|20bjl and (|36| and neglecting the decay- 
ing mode, we can write the exact nonlinear growth fac- 
tor for the density inhomogeneity in our model. Defining 
D + = 6 /Si, we find 

n z l f + {l-5^z l B[{x + 1 f+ [y + uf]} 
+ l-S9 Zi {f + + B[(x + 7 ) 2 + (y + ujy]y [ ' 



which, for oj = r ) = x = y = Q, simplifies to 



z%f-t 



5° 1-SUiU 



(52) 



where zi is the redshift of the initial condition Si and B is 
defined in Eq. (|23|) . We have substituted the dependence 
on f3 + for — Sfz.i, where Sf is the initial density pertur- 
bation along the z-axis. In our model we are only free 
to choose the distribution of S along the z-axis, with the 
distribution of the density along the x- and y-axis then 
being given by the metric. This is easily understood by 
considering the fact that all free functions in the metric 
are only function of the z-coordinate. The only function 
in the metric containing x and y is A with its dependence 
being fixed (see Eq. ((30]) ). 



section are the expansion scalar, the shear tensor and 
the electric part of the Weyl tensor. They are defined as, 
respectively, 



i 



Vab = U(o;6) ~ 7^>h a b + U(a u b)i (53b) 



F a b — C ac bdU C U d , 

where u a is the fluid 4 velocity 

u a = [-1,0,0,0], 

and we define 

Ua = U a ;bU b , 
hab = 9ab + U a U b 



(53a) 
(53b) 
(53c) 

(54) 

(55) 
(56) 



and Cabcd is the Weyl tensor. 

From the above definition and using the metric (|15D . 
we find 



6 = 3- 



Z 



(57) 



To analyze the shear and Weyl tensors, we need some 
definitions. As expected [49| . from the metric we find 
that H a b = ^>ab = and hence we are considering mod- 
els in the "silent universes" class [lH, i.e. cosmological 
models where there is no communication between fluid 
elements. Since the equations evolving p, 0, <j a b and E a u 
are ODEs, each fluid element evolves independentl-yT"!. 
For an analysis of the dynamical systems for these co- 
variant variables in "silent models" , see [46| and [47 1, for 
the case with and without A, respectively (see also 48j). 

Exact solutions in the "silent universes" class include 
Bianchi I and Szekeres models, as proved in [HI]. The ex- 
istence of more general "silent" exact solutions of Petrov 
type I is unlikely for A = (see [731 and references 
therein), but few exist for A ^ [74[. However, the 
"silent" approximation with H ab = u ab = holds true 
for first-order scalar perturbations [65| and at second or- 
der outside the horizon (75[. It also corresponds, in the 
covariant description, to the so called long wavelength 
approximation in gravitational collapse (see [76| and ref- 
erences therein), with H ab only becoming non- negligible 
in between different Kasner phases l53l|. 

In "silent universes" , according to 49i] , a ab and E a b 
have a common cigenframe. Wc can thus expand them 
as 



V. COVARIANT VARIABLES 

A. Deriving the variables from the metric 

We now consider the covariant fluid flow description of 
our space-times [M US S IM 0] • For the metric (JTSJ) 
the magnetic part of the Weyl tensor H ab is known to 
be zero [49| and our dust flow is irrotational, u! ab = 0, 
therefore the only variables wc need to consider in this 



E a b — ^ ' E a e aa G a b, &ab — ^ &a£aa£ab- (58) 



a=l 



a=l 



The orthonormal tetrad e aa can be found to be e\ a = 
SSl, e-ia = SS\ and e 3a = SZS\. From this tetrad, we 



10 Obviously, initial conditions have to satisfy spatial constraints, 
but here we are only concerned with the local time evolution. 
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can now find the complex null tetrad required to compute 
the Weyl scalars in Section BlI CI The procedure outlined 
in [77} has been followed and the complex null tetrad we 
find is 



S S 
m a = (0, -j=,-i-j=,Q), 

1 SZ 
1 S 1 / 7 

*° = ( ~7!' ' ' 



(59a) 
(59b) 
(59c) 



and rr% a is the complex conjugate of m a . Returning to 
the shear and electric part of the Weyl tensor, they are 
trace free and so we write 



(60) 



Q = l 



We can then define a + = \{<J\ + 02), <J- = -^=(a\ — 
<r a ) and E+ = \{E X + E 2 ), E_ = ^{E 1 - E 2 ) as a 
convenient set of dynamical variables, which determine 
the dynamics of the shear tensor and the electric part of 
the Weyl tensor completely. Using these definitions, we 
find 



and 



E+ = 



1 Z 

~3~z' 

M 



0, 



(61) 



Po _ 
6S 3 6S 3 Z 



1 



= -qPS, £-=0, (62) 



where in the last equation for E + we also used the EFEs. 
Comparing with Eq. (|43|) we see that ^ = — E+. 



B. Phase plane analysis 

The complete set of ODEs governing the dynamics of 
our models is given by 



at - E. 



+ • 



E+ = -QE + - 3a + E+ - -pa, 
P = -®P, 

e = -ie 2 -64-i P + A. 



(63a) 

(63b) 
(63c) 
(63d) 



This system of ODEs is the subset, for cr_ = = 0, of 
that considered in [47], 13 f° r "silent models" with A. 

As in the previous sections, our aim is to decouple the 
dynamics of the background from that of the inhomo- 
geneities. The problem is that the variables p and 
incorporate both a background and an inhomogeneous 
part. On the other hand, the shear and the electric Weyl 
tensor vanish in an FLRW space-time; for this reason 
they are first-order gauge-invariant variables HD, |7S| 



and here represent exact inhomogencities. For O and p 
we can write 

O = 6 + 6, (64) 
p = p(l + S), (65) 
where, using Eq. (|57[) . we define 



and 



9 = s| = 3H, 



Z 

z' 



(66) 



(67) 



and from comparing with Eq. ()61|) . we find 8 = — 3tr+. 

Given expressions for 9 and a+ from the line element, 
one can use Eqs. (|63p to find the evolution equations 
for S and Z . Also, in finding E + there is an integration 
constant, which corresponds to po, yet here one finds this 
function to be space dependent in general. However, we 
can find another expression for E + from the differential 
equation for a + . We find 

2 SZ IZ 
E + = 3SZ + 3Z- ^ 

The aim here is to describe the full dynamics of our 
models with as few variables and differential equations as 
possible. We will now need to split the covariant variables 
into background and inhomogencities. To analyze the 
system, we start from the Friedmann equations for a flat 
universe; this will eventually reduce our system by one 
order. We start from 

H + H' 2 + -p - \k = 0, 



6 r 
H 2 



1 



1 



-A. 



(69a) 
(69b) 



Combining the two equations, one finds the continuity 
equation for the background (see also Appendix IA"|) 



P 



-3Hp. 



Changing variables to 



0/ 



A 
3IP' 
P 

3H 2 ' 



we can rewrite Eq. (|69b[) as 

1 = Q m + Q\. 
Then, using Eq. ()69a|) . we find 

n' A = 3n A (i - n A ) , 



(70) 

(71a) 
(71b) 

(72) 
(73) 



where ()' has been defined as ()' = -g°(). Eq. ([73]) encodes 
the dynamics of the background variables. From the def- 
initions (|64|) and (|65|) and continuity equation (|63cp . we 
obtain 



5' 



H 



(1 + 6). 



(74) 
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Rccxpressing the evolution equation for er + (|63a[) in the 
new time and background variables we obtain 



3 2 
r + 

We now define the new variable 



-2 -+ + ^ + y(l-n A )& (75) 



(76) 



which represents an expansion normalized velocity de- 
viation. Note here that this variable S + is different to 
the one used in reference (48[, where the normalization 
of cr + is given by 0. We have chosen this normaliza- 
tion as it leads to a separation of background and per- 
turbation variables. In addition, H being a monotonic 
function describing the expansion of the background, the 
time derivative ()' = is uniquely defined across turn- 

I 



around, when 9 = — O. Now, consider the system 
Q' A = 3f2 A (l-0 A ), 

8' = 3E+(l + a), (77) 
£ + = i£ + (3ft A + l)-3E2_-i(l-fi A )<5. 

We can further reduce the order of this dynamical system 
using the conserved quantity d = j3+B: using Eq. f) 13[) . 
we can find the constraint 



1(0 A -1)<5-|±(5 + 1)(^ 



/2 _ Q 3/2 
A "A 



2/3 



. (78) 



where j3+ and A have been defined in Eqs. (fl4j) and (p~0|> 
respectively. Using this constraint, we can decouple the 
two differential equations for 8 and £+. Then, our dy- 
namical system takes the final form 



a =30 A (i-n A ), 



8' = - 



(i-n A )s+e±(8- 



,3/2 



2/3 



(1 + 5). 



(79) 



Therefore, we have reduced the system of four coupled 
differential equations to a system of two differential equa- 
tions, where the background evolution has decoupled 
from the deviations. The reduction in order has been 
due to the assumption that K = and the result that 
the curvature deviation is constant in our space-time, 
d = 0. The decoupling has been made possible by choos- 
ing a suitable set of variables to work with. The system 
is now fully characterized by the differential equations 
for fl A and 8. £+ can then easily be found using the 
constraint equation (|75j) . The dynamics of the variables 
can be shown in compressed form by doing a phase plane 
analysis. Plots of this analysis are shown in Fig. ([4]) for 
the simple case of undcr-densities. In the plots we have 
suppressed the decaying mode of F for clarity. 

In both the figures, the initial conditions were set very 



close to the 0/ 



E + = point and again for 



clarity, we have plotted a single trajectory for each value 
of /3+. This point is associated with Einstein-de-Sitter 
(EdS) space. The different evolutions are due to different 
values of /3+ and all trajectories in both figures tend to de 
Sitter (dS) space. In the first plot, dS space is represented 
by all values of 8 along the fi A = 1 line, i.e. each point on 
this line is a fixed pointF^I. In the second plot dS space is 
just represented by the fl A = 1, £+ =0 point. Therefore 



11 Notice that when <5 is a constant, the total density p scales with 
the background p ~ S~ 3 and therefore tends to zero at late 
times, cf. |44H for a general second order analysis of the asymp- 
totic evolution of perturbations in ACDM and the cosmic no- hair 
theorem that follows. 



we see that dS is always the late time attractor for all 
solutions, consistently with the analysis in [43| . 



VI. CONCLUSION 

The cosmological constant problem has lead to the 
investigation of alternatives to the concordance ACDM 
model, considering other forms of dark energy or the- 
ories of gravity alternative to GR. However, the con- 
sequences of Einstein's equations and their nonlincarity 
have not yet been fully explored in cosmology. We typ- 
ically use a homogeneous isotropic FLRW background 
plus perturbations, treating nonlinear structure forma- 
tion with the Newtonian approximation. Most observa- 
tions are interpreted assuming this Friedmannian frame- 
work, even completely neglecting inhomogeneities, as we 
normally do in determining distances on the base of a 
FLRW distance-redshift relation, cf. [Tl |. 

So far there is good observational support for a ACDM 
cosmology build in this way, when compared with al- 
ternative theories of gravity under the same Fried- 
mannian assumptions (sec e.g. [7^81[), i.e. using a 
post-Friedmannian approach (see 82h85| and references 
therein). At a time when galaxy surveys and other ob- 
servations are reaching unprecedented sky coverage and 
precision it seems however timely to investigate, in par- 
allel with this linear cosmology approach to alternative 
gravity theories, the effects of nonlinear general relativis- 
tic dynamics on the growth of structures and on observa- 
tions. Analytic inhomogeneous cosmological models are 
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0.15 



FIG. 4. Phase planes for the inhomogeneity variables 5 and 
E+ versus the background variable JIa, for the case of under- 
densities, /3+ > 0. In each plot we only consider the growing 
mode and, for clarity, we plot a single trajectory for each value 
of/?+. 



indispensable to analyze and understand these effects in 
a simplified context. 

Assuming GR, in this paper we have found and ana- 
lyzed exact solutions of Einstein's equations describing 
an inhomogeneous universe with pressureless dust (rep- 
resenting CDM) and a cosmological constant A. Our 
models are of the Szekeres class II type, generalized to 
the A ^ case; these solutions were previously obtained 
in the original Szekeres form of the metric by Barrow 
and Stein-Schabes [43j |. We obtained our models us- 
ing the metric form previously used for the A = 
case by Goode and Wainwright j45(. The great advan- 
tage of this form of the metric is that it allows to split 
the dynamics of the model into a part that describes a 
ACDM FLRW background, which we have taken to be 
flat, and a part describing an inhomogeneous deviation 



from this background, represented by a single function 
Z in the line element (p}. The coordinates in the metric 
(|T|) are synchronous and comoving, the relativistic analog 
of Newtonian Lagrangian coordinates, thus Z represents 
the extra change - with respect to the general Hubble ex- 
pansion - of the physical distance between pairs of fluid 
elements along the z-axis. Einstein's equations dictate 
that Z = A(x) + F(t, z), where the dependence of Z (A) 
on the x, y coordinates is fixed, while the dependence on z 
is arbitrary. The density deviation from the background 
has the simple form 6 = —F/Z. This can be then used to 
model an arbitrary initial matter distribution along one 
line of sight. The coordinates in the metric (QJ are syn- 
chronous and comoving, the relativistic analog of Newto- 
nian Lagrangian coordinates, thus Z represents the extra 
change - with respect to the general Hubble expansion - 
of the physical distance between pairs of fluid elements 
along the z-axis. 

A very useful property of the Z function is that its time 
dependent part, F, satisfies the same linear second order 
ordinary differential equation (|12jl that 6 satisfies at first 
perturbative order in Newtonian theory Q. This same 
equation is satisfied by a A variable in gauge-invariant 
perturbation theory [64l468l | , where A reduces to S in any 
comoving synchronous gauge (see (oTI . l68j ) . This linear- 
ity property of Z implies that it satisfies a superposition 
principle, i.e. we can build an arbitrary initial matter dis- 
tribution along the z-axis, e.g. adding up Fourier modes, 
and define a spatial average. The Z inhomogeneity does 
not need to have zero average along z, thus in this sense 
the FLRW background may or may not be representa- 
tive of an average. This could be used to study a simple 
example of averaging leading the a redefinition of the 
background; we leave this and other issues for a future 
analysis [4J|. For instance, also worth of study is that 
the linear superpositions of modes in Z leads to mode 
coupling in the nonlinear d — —F/Z. 

The two linearly independent solutions for F are the 
well known growing and decaying modes for <5 in pertur- 
bation theory @]. We have shown that a crucial quan- 
tity is a conserved curvature inhomogeneity, which is a 
first integral of the dynamics: this generates the growing 
mode, in complete analogy with perturbation theory. As- 
suming a vanishing decaying mode, we have given an ex- 
act formula for the nonlinear growth factor in our model. 

We have also studied the local formation of singulari- 
ties in our model, finding in particular that, at any given 
comoving point (i.e. for a fixed fluid element), for S > 
pancake may form, similar to the Zel'dovich pancakes of 
Newtonian theory. Unlike for the A = case, pancakes 
are not unavoidable, given that the growing mode in Z 
for A 7^ tends to a constant value. Therefore, at any 
given point in our synchronous comoving coordinates, ini- 
tial conditions may be found such that the pancake never 
forms. 

The models studied here belong to the "silent" class 
HI with A studied in (cf. [H and (H). A patch in 
our model not evolving to a pancake unavoidably asymp- 
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totically approaches de Sitter [46| (cf. [^^^satisfying the 
cosmic no-hair theorem. With respect to [U EH we have 
greatly simplified the dynamical system analysis, intro- 
ducing new variables that lead to a decoupling of the sys- 
tem, so that two variables only are needed: describing 
the background and 6 (or a shear variable) describing the 
inhomogeneity. 

It is known [49| that the Szekeres metric is of Petrov 
type D: since this is an algebraic property of the geome- 
try (the metric), it is unaffected by the particular form of 
the field equations. We have computed explicitly ^2, the 
only non- vanishing Weyl scalar (in a specific null tetrad) 
for Petrov type D space-times, showing that it does not 
depend on A. The issue of the contribution of A to lens- 
ing has been a subject of debate recently, see [63| and 
references therein. \J> 2 is the only Weyl contribution to 
the geodesic deviation equation and A contributes to the 
FLRW background part of the metric but not to the in- 
homogeneous part of the Ricci tensor. It then follows 
that A can only contribute to lensing indirectly, through 
its effect on the background expansion. 

Finally, it is worth noticing that in the setting we use, 
of a FLRW background plus exact inhomogeneous devia- 
tions, the question arises of the gauge-invariance of these 
deviations, and the relation that exist to the perturbative 
gauge-invariance (or variance) . Clearly, we have used a S 
that is defined in our synchronous comoving coordinates: 
while p itself is a scalar field on the inhomogeneous space- 
time and as such is invariant, defining d requires a map 
between the background and the inho mog eneous space- 
time. As it was shown in [8(| (see also |87|). it turns out 
that scalars that vanish in the background, and that are 
therefore gauge-invariant at first-order 164 l88l . HH, are 
also gauge-invariants of the exact theorjo- The Weyl 
scalar ^2 mentioned above is therefore an exact gauge- 
invariant description of the deviation between the back- 
ground FLRW and the exact inhomogeneous space-time. 

Having developed and analyzed our model, we would 
like to use it in a future work [4l| to study more in the 
detail the growth of nonlinear density inhomogeneities 
and light propagation through them, to see how obser- 
vations may be affected. In particular, we would like to 
investigate whether the distance-redshift relation is the 
same in an inhomogeneous space-time as in FLRW or 
whether non-negligible corrections need to be made. In 
this respect, the aim will be to produce an analysis of 
this problem based on exact solutions with a continuous 
inhomogeneous distribution of matter; such a study will 
be complementary to the work by Clifton and Ferreira 
[HJ, who used an approximate model where all matter 



is concentrated in a pointlikc lattice distribution. Our 
modcls arc certainly well suited for this task, as we have 
shown that nonlinear structures can be modeled along 
one line of sight without having to resort to any approx- 
imations. Clearly, our model can also be used to study 
the growth of nonlinear structures. Indeed, at least at 
an initial time, we can specify an arbitrary matter dis- 
tribution along one direction. Simulating the growth of 
large over-densities can result in pancake singularities in 
the density field after a certain time. This is not an issue 
though, if one tries to simulate realistic density profiles, 
as we find that - in a typical ACDM cosmology - this 
problem only occurs when the perturbations during the 
epoch of recombination are set to be much larger than 
10~ 4 . Therefore simulating structure formation and per- 
forming light tracing along one line of sight using realistic 
initial conditions, based on a typical perturbation power 
spectrum, is well within the scope of our model. 



Appendix A: The continuity equation 

In this Appendix we look at the continuity equation, 
at first with no assumptions on the theory of gravity and 
the field equations, but assuming the line element (fTJ). 
First of all, let us derive the general form of p. We must 
solve 



-Qp, 



(Al) 



where an expression for in terms of metric functions is 
given in Eq. (|57p . Using this expression, we modify this 
equation to find 



(S 3 Z) 
S 3 Z 



(A2) 



Hence we can find the general solution for the density 

M(x) 



P = 



S 3 Z 



(A3) 



which is exactly the form stated in Eq. 

Given that S = S(t), we now assume that p can be 
written as 



p = p(l + S), 



(A4) 



where p = p{i) is assumed to be the homogeneous density 
of a FLRW space-time with scale factor S(t). Under 
this assumption p = p~oS~ 3 and 8 is the dimensionlcss 
density deviation from the background p. Substituting 
this decomposition of p into the continuity equation (|Aip 
yields 



12 An important subtlety is that in general a field that vanishes in 
the background, and that is thus gauge-invariant at first-order, is 
not gauge-invariant at higher order, although is gauge-invariant 
for the exact theory: perturbations "live in the background" , the 
exact deviations are fields on the inhomogeneous space-time, see 
~ and M 




= 0, 



(A5) 



where we identify p = — 3§p as the background conti- 
nuity equation. Subtracting this off and rearranging, we 
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obtain 



which we used in Eq. (|34l) . 



1 + 6 



Z 

z 



(A6) 



Appendix B: Details on solving the EFEs 



Integrating this equation we find 

. (1 + Si)Zi - z 



(A7) 



where Si and Zj are arbitrary initial values. 

If we now assume the EFEs, it follows that Z(t, x) = 
F(z, i)+A(x), sec Eq. ([B3|> below, such that Zi =F t +A. 
With this, we can substitute into (IA7I) to obtain 



In this Appendix we solve EFEs ([2]) for the line ele- 
ment ((T|) with the energy momentum tensor given in 
Since we are only considering irrotational dust and a cos- 
mological constant, we can use the off-diagonal terms of 
G a b as constraints. We start by finding 



G, 



G t y — 



z 



z 



1)1 



0. 

0. 



(Bl) 
(B2) 



F + A' 



(A8) 



which is exactly the same form of 6 as presented in Eq. 
(l36|) . Note that in obtaining (|A8|) . we have assumed 



Si = - 



F, 



Fi + A 



(A9) 



Knowing that F has a decaying mode /_ = /_(t) and 
a growing mode / + = f+(t) (see Appendix [B] below), 
with the latter such that /+(0) = 0, we have chosen the 
relation between initial conditions (|A9[) such that, in the 
case of a purely growing mode, (5(0) = 0. Vice versa, 
assuming 6(0) = implies Eq. (|A9j) . 

Finally, with the choice §A9§, it follows from (|A8jl . (|A3]l 
and (TAHl that 



M = p A 



(A10) 



This implies that we can write 

Z(x,t) = F(z,t)+A(x) 
Furthermore we find 

2a z Z t 



= 0, 



(B3) 



(B4) 



and since we know that Z will change with time we find 
that 

a(x) = a(x, y). (B5) 
The last restriction from off-diagonal terms is given by 
the G xy term. It dictates 



a y Z x + a x Z y 



0; 



(B6) 



we will use this constraint later. The four diagonal terms 
give the equations 



Gtt — 

Gxx _ n S t t 
S 2 - S 

Gyy _ Su 
S 2 s 

G zz 



O St Zt i Q f St' 

z s z +t) \s, 



S 2 Z 2 



S 2 e 2 



'St' 



o St_ Zt_ I Ztt 

° s z ~ z 



o S t Z t 

s z 



Zt, 

z 



{fit + ^r + a xx + «„„) = p + A, 
^~ ZS 2 e 2a ^P^y^v Ot-x^tx ^yy) — ^-j 
\~ zS^e la ( a x^x ~ WyZy ~ %xx) = A, 



2^ + (f) 



e 2a (p^yy Q-xx) 



= A. 



(B7a) 
(B7b) 
(B7c) 
(B7d) 



r 



From this system, one can see that by subtracting Eqs. 
(|B7b[) and (|B7cj) from each other, one obtains 



2a x Z x 



2a y Z y 



Z„ 



0. 



Multiplying Eq. (|B7d|) by S 2 S t gives 



2SttSSt + S t 



-(a xx + Oiyy) = AS 2 S t 



(B8) 



(B9) 



time. Integrating this equation with respect to time then 
gives 

SS?--^(a xx + a yy ) = ±AS 3 +C, (B10) 

where C is an integration constant. Rearranging this 
equation one finds 



Together with Eq. (|B5[) . given that S = S(i), this shows 



that the term e 2a (a xx 



■a v 



is a constant in space and 



S t \ C 1 . K 



SJ S 3 + 3 A+ S 2 ' 



(Bll) 
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which is the Friedmann constraint equation for ACDM. 
We have defined 



K 



(B12) 



the constant term found in (|B9[) . This differential equa- 
tion admits the solution 



1 



l + ±if(x 2 + y 2 )' 



(B13) 



Thus we find for K = that e a = 1 and therefore a x = 
a s = = <x yl/ = 0. From now on we restrict our 

attention to a flat universe. 

In the Friedmann equation (|B11[) . we identify the term 
C/S 3 with the energy-density term p/3. This satisfies the 
background continuity equation (see Appendix [A} that 
gives 

Po 



(B14) 

and therefore we identity our constant C as 

3C = p . (B15) 

Since we are only considering a flat universe, we find 

Z xy = 0, (B16) 

and 

Z X X = Zyy, (B17) 



from Eqs. (|B6|) and (|B8|) respectively. We can thus reduce 
the system Eqs. (|B7a) - (|B7d)) to 



2 5 Z +i 



- 2 



ZS 2 



= P 



s 



4+ 



Zyy 

zs 2 

St 

s 



A, (B18a) 
A,(B18b) 
A. (B18c) 



The combination (|B18b|) + 4 (|B18c|l - jr (|B18a|) gives 



Clearly the LHS is a function of time and z only and the 
RHS is a function of x, y and z only. Hence both sides 
must be equal to a function of z only. Call this function 
g(z): then F satisfies the differential equation 



2§« 



Po 
2S* 3 



F 



9M 

S 3 ■ 



(B23) 



This ODE has two homogeneous and one particular so- 
lution. We denote the homogeneous solution F = 
F h (z,t). The particular solution is easily spotted to be 



F" 



2g(z 
Po 



This gives the form 



F = F h (z,t)- 
Also, from Eq. (|B22l we find 



.4 



2g(*) 
Po 



2g{z) 
Po 

M 
Po ' 



which implies 



F h (z,t) + 



M 
Po ' 



(B24) 



(B25) 



(B26) 



(B27) 



This shows that the metric is completely independent 
of the function g(z). This can be understood by looking 
back at Eq. (|B3|) . There we have decomposed the func- 
tion Z into two separate functions A and F, but both 
those functions are functions of z, which means that there 
is always a certain arbitrariness in the choice of F and 
A. We could easily add g(z) to F and subtract it from A 
and still end up with the same function Z . Thus, we can 
choose here the function g{z) to be equal to zero without 
loss of generality, which gives 



A=—. 

Po 



(B28) 



Stt Ztt 2 StZt__P A 
6 S + Z +2 S Z - 2 +A ' 

Combining Eqs. QBlTJ) , (|BT5| and (|B18cj) yields 



Stt 

~s~ 



Po 
6S 3 



(B19) 



(B20) 



which is the Friedmann equation for ACDM. Subtracting 
3 times (|B20j) from (|B19[) gives 



Ztt r-St Zf 

~z~ ~s~z 



M 
2S 3 Z 



«L =0 
25 3 ' 



(B21) 



where we have used Eq. (|A3|) . Substituting the decom- 
position of Z we found in Eq. (|B3[) . we find 



S 3 F ti 



2S 2 S t F t 



Po 



F 



M 

Y 



Po 



(B22) 



Keeping the function g(z) would not change any results 
but would merely clutter the equations. In other words, 
we only need the homogeneous part of Eq. (|B23|) to com- 
pletely specify our solution. From now on therefore, we 
will assume that F satisfies the equation 



Ft, 



2^F t 
S 



Po 
2S 3 



F = 0. 



(B29) 



As first pointed out in [45[, this linear ODE is the equa- 
tion satisfied by first-order density fluctuations in a per- 
turbed dust FLRW universe. 

Since this equation has two linearly independent solu- 
tions, which exhibit growing and decaying behavior, we 
can write F as 



F(z,t) = P+(z)f+(t) + P-(z)f-(t), 



(B30) 
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where /+ (t) and /_ (t) are the growing and decaying solu- 
tions respectively and (3+(z) and ft-(z) are free functions 
of z. Using this decomposition of F and Eqs. (|B16|) and 
(|B17[) we can find the functional form of A, obtaining 

A(x) = a(z) + b(z)x + c{z)y + d{z)(x 2 + y 2 ). (B31) 

Looking at the system of differential equations, Eqs. 
(|B18ap ~ ()B18c[) . we can see that we have only extracted 
two equations out of this system, so we should be able to 
get more information out of it. Subtracting Eq. (|B18bl) 
from Eq. (|B18c|) we find 



„St Z t Zu Z xx _ 
~S~Z + ~Z~~ZS 2 ~ ~ ' 



(B32) 



which can be brought into a slightly different form by 
using the relation between Z and F, 



3^F t + F a 



S 2 



0. 



(B33) 



Using Eq. (|B29[) and noting that Z xx = 2d, we find 



(B34) 



This equation really takes the form of a first integral 
equation of Eq. (|B29|) . see Appendix [Cl below. We discuss 
the form of this equation in more detail in Sec. IIV1 on 
perturbation theory. 



Appendix C: The relation between the first and 
second order equations for F 

We have found two different differential equations for 
F, one first-order (|B34[) and one second order (|B29|) . We 
explicitly show here that the first-order equation is the 
first integral of the second order equation when K = 0. 
We start by modifying Eq. (|B29[) to obtain 



St 

s 



Ft, 



Fi 



S 2 



F = 0. 



From Eqs. flBlTJ) and (|B20|) we can find 

2 - ci / ci \ 2 



3 

s 



p 

2 



Stt 

s 



(Cl) 



(C2) 



Using this expression, we can rewrite Eq. (|C1[) as 
SStFtt + SSttFt + StStFt 



2S~ " 2S 2 
which can easily be integrated to obtain 



P ° Ft - ^FS t = 0, (C3) 



S 1 2S 3 S 2 



0. 



(C4) 



where C is a constant in time. This equation takes the 
same form as Eq. (|B34j) . However the function C here 
cannot be related to the metric, whereas, when this dif- 
ferential equation is derived using the EFEs, we can find 
that C = d(z), which features in the metric itself. 
Appendix D: Dimensional Analysis 

We add this Appendix to aid the reader in gaining 
some physical interpretations of some of the variables. 



In this section L n denotes the dimension of length to the 
power of n; assuming c = 87rG = 1, we find 



[Z] = L°, 


(Dl) 


[*] = L, 


(D2) 


U]=L-\ 


(D3) 


x \ = [y] = \A = L , 


(D4) 


IS] = l°, 


(D5) 


\n 1 — T~ 2 




[A] = L- 2 , 


(D7) 


[B] = L- 2 , 


(D8) 


[r]=L°. 


(D9) 



Appendix E: Symmetries 

We would like to show here that for certain choices of 
the free functions of the metric, we find axial symmetry. 
From the line element Eq. (|15[) . we find the metric 



9ab 



( -1 \ 

S(t) 2 

S(t) 2 

V S{t) 2 Z{t,x,y lZ ) 2 J 



(El) 



Clearly this metric is dependent on all 4 space-time vari- 
ables. We can decompose Z(t, x, y, z) as 



Z = 1 + F(z, t) + BP+ [{x + 7 ) 2 + (y + uif 



(E2) 



where f3+ = f3+(z), 7 = 7(2) and to = lo(z). If we choose 
7 = tjj = (choosing 7 = Ci and uj = C2 gives the same 
result) we find 

Z{x, y,z) = l + F(z, t) + B(3 + (z) (x 2 + y 2 ) 



By making the coordinate transformation 
x = psin(</>), y = pcos(^), 
we can rewrite Z as 

Z(p,z) = l + F(z,t) + B(3 + (z)p 2 . 
In the new coordinates (t, p, (j), z) we find 

/ -1 

S(t) 2 
S{t) 2 p 2 



g a b 



V 







(E3) 
(E4) 
(E5) 

(E6) 



S(t) 2 Z(t,p,z) 2 J 



The metric does not depend on <p, hence, for 7 = ui = 
the solution has an axial symmetry about the z-axis. 
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